Generalized contact process on random environments 
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Spreading from a seed is studied by Monte Carlo simulation on a square lattice with two types of sites 
affecting the rates of birth and death. These systems exhibit a critical transition between survival and 
extinction. For time-dependent background, this transition is equivalent to those found in homogeneous 
systems (i.e. to directed percolation). For frozen backgrounds, the appearance of Griffiths phase prevents the 
accurate analysis of this transition. For long times in the subcritical region, spreading remains localized in 
compact (rather than ramified) patches, and the average number of occupied sites increases logarithmically 
in the surviving trials. 

PACS numbers: 02.50.-r, 05.50.+q, 87.23. Cc 



I. INTRODUCTION 

Contact process (CP) was introduced by Harris jj| to 
describe competition between death and birth events for 
a spatially distributed species. In this model, the or- 
ganisms localized at the sites of a lattice can die with 
a death rate or can create additional offspring in one of 
their (empty) neighboring sites. When increasing the ra- 
tio of birth to death rates, this system exhibits a critical 
phase transition from empty to the active state as the sys- 
tem size approaches infinity. This transition belongs to 
the so-called directed percolation (DP) universality class 
[^|,^). Similar transitions have been observed for some 
other, related models (for recent reviews see the Refs. 
ID- 

The above mentioned CPs applied homogeneous back- 
grounds. Numerical investigations on two-dimensional 
random environments have been restricted to diluted 
lattices p|-|lO|]. Noest have shown that, although 
the transition remains continuous, the extinction process 
is modified drastically on the randomly diluted lattice 
above the percolation threshold. Recently Dickman and 
Moreira || have reported that the scaling (a fundamen- 
tal feature of critical transitions) is violated in this case. 
According to a field-theoretic analysis, similar result is 
predicted by Janssen JTT[ ] . 

Below the percolation threshold on a diluted lattice, 
the available sites form isolated (finite) clusters. In any 
particular cluster the species is expected to die out within 
finite time and the cluster remains empty afterwards. A 
simple mathematical description suggested by Noest |?]] 
indicates that the average concentration of species van- 
ishes algebrically with time and this phenomenon is anal- 
ogous to the relaxation observed for the "Griffiths phase" 
in disordered spin models 0] . 

Field-theoretic arguments [0 support that the tem- 
porally quenched disorder itself, as well as the spatially 
quenched disorder, disturbs crucially the DP transition. 



However, DP transition is expected on large scales for 
those random environments where the spatio-temporal 
disorder is uncorrelated. For two-dimensional systems 
this expectation is supported only by a small number 
of numerical evidences. For example, in an evolution- 
ary game, extinction of one of three competing strategies 
represents a DP transition on a background whose time 
dependence is maintained by the competition of the re- 
maining two strategies Jlj|. 

In the present paper, we study a generalized CP on a 
square lattice with two types of sites that provide differ- 
ent conditions for survival and reproduction for a hypo- 
thetical species. Our analysis is not restricted to frozen 
backgrounds. We give numerical evidence that the ex- 
tinction process becomes analogous to DP transition if 
the random environment involves uncorrelated time de- 
pendence. Using numerical simulations, we investigate 
the main characteristics of spreading when initially there 
is only a single individual in the system. This method 
was suggested by Grassberger and de la Torre jl4| to 
study DP transition (for homogeneous background) in 
the close vicinity of the critical point. Moreira and Dick- 
man H have demonstrated that this technique is also 
efficient for the investigation of inhomogeneous systems. 
Now we have adopted this method to study the above- 
mentioned generalized CP in random environments. The 
results confirm the theoretical expectations |^,|^|ll],[l5| as 
well as the previous observations based on MC simu- 
lations J6|— |l0|] . We have also studied some geometrical 
features of the patches formed by occupied sites in the 
subcritical region. 



II. THE MODEL 

Space is represented by a square lattice. Each lattice 
site, r = (x, y) {x and y are integers), is assumed to repre- 
sent a microhabit for one individual of a species. Quality 
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of a site, g(r), determines the probabilities of survival and 
reproduction of the inhabitant individual. We assume 
two microhabitat types: good and bad (where g(r) = 1 or 
0, respectively). Site qualities are distributed randomly 
in space. Initially, each site is chosen to be good (bad) 
with a probability P (1 — P). Distribution of individuals 
over this random background is described by a state vari- 
able <r(r) that is 1 for occupied and for empty state. 
Time evolution of the system is governed by three ele- 
mentary processes affecting the values of cr(r) and g(r). 
For an occupied site, extinction (death) can occur with a 
death rate c%(r)]. For an empty site, colonization (birth) 
can take place with a birth rate 6[<?(r)], provided that a 
randomly chosen nearest-neighbor site is occupied. It is 
supposed that good sites provide better conditions for 
living, that is 6(1) > 6(0) and d(l) < d(Q). Furthermore, 
redrawing the value of g(r) we change the environment 
with a rate /. The new state is independent of the actual 
values of g(r) and <r(r). This occasional modification re- 
sults in stationary distribution of good or bad sites over 
time, with probabilities P and 1 — P, respectively, and 
leaves their distribution uncorrelated. 

In the special case when every site is good (P = 1), the 
present model is equivalent to the well investigated con- 
tact process on homogeneous background . Previous 
studies have shown that in this case, the species dies out 
if A = 6(l)/d(l) > A c = 1.6488(1). Concentration of 
population vanishes algebrically in the active stationary 
state, that is c oc (A - \ C Y where (3 = 0.575(3) if A -> A c 
P,^| Jl^Jl7|] . This critical transition is accompaned with 
diverging fluctuations and correlation length. 

The conditions of CP on diluted lattice can be repro- 
duced by allowing the species to stay only in good sites 
(6(0) = and d(0) = oc) for P < 1 and / = 0. 

A mean-field approximation for frozen backgrounds 
(/ = 0) predicts different concentration of species on 
good and bad sites: c\ and cq, respectively. Neglecting 
the details of a straightforward calculation, the concen- 
trations obey the following expressions: 
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(1) 
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where X a — b(a)/d(a) (a — or 1), and the average 
concentration, c = Pc\ + (1 — P)cq is obtained as 



A+ y / A 2 + 4A A 1 (A av -l) 
AqAi 



(2) 



where A = A Ai — A — Ai and A av = PAi + (1 — P)Xq- 
This mean-field approximation predicts that concentra- 
tion vanishes linearly when A av — > +1. More precisely, 
for low concentrations the above expression can be ap- 
proximated as 



PAi + (1 - P)A - 1 
Ao + Ai — AqAi 



(3) 



if A av > 1. Obviously, the trivial (c = 0) solution will be 
stable for A av < 1. 



In the limit / — > oo death and birth occur with aver- 
aged rates: d = Pd(l) + (1 - P)d(0) and 6 = P6(l) + (1 - 
P)6(0)), respectively. Thus the system becomes equiv- 
alent to a homogeneous system (where Co = c\ = c), 
and the average concentration obeys the traditional form 



A- 1 
A 



(4) 



if A = b/d > 1 and c = otherwise. The difference 
between the Eqs. (0) and (ji|) refers to the / dependence. 

III. SIMULATION OF SPREADING 

The Monte Carlo (MC) simulations are performed on 
a lattice with L x L sites under periodic boundary con- 
ditions. Followingthe method suggested by Grassberger 
and de la Torre Q spreading is investigated by averag- 
ing over many trials M when the initial state is close to 
the absorbing state. More precisely, each run is started 
with a single individual [at position r = (0, 0)] on a new, 
random, uncorrelated background. We have determined 
the survival probability 



S(f) = <*(»(*))> > 



(5) 



where the number of individuals at a given time t is de- 
fined as 



i(t) = £>(r) 



(0) 



and (. . .} means the average over M trials and 8(z) — 1 
(0) for z > (z < 0). We have also evaluated the average 
number of surviving individuals, 



N(t) = (n(t)> 



(7) 



and the mean-square distance of individuals from the ori- 
gin, 



R 2 (t) = 



l 



N(t) 



^r 2 a(r) 



(8) 



M is varied from 10 4 to 10 8 in the systematic inves- 
tigations. The system size L is adjusted to exceed sig- 
nificantly the average radius of occupied sites, namely, 
L <^ l5R(t m ax), where t max indicates the time limit. As 
usual, time is measured in Monte Carlo Steps (MCS). 
During one time unit, each site has a chance to modify 
the value of tr(r) (approximately once on average). 

The efficiency of these simulations can be greatly im- 
proved by labeling the individuals, recording their co- 
ordinates, and picking up one individual randomly for 
updating the state of occupancy || . In the present case, 
a chosen individual residing at site r dies with a prob- 
ability d(g(r))/(b(l) + d(0)) or creates an offspring on 
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one of its nearest-neighboring sites r' with a probability 
b(g(r')) / (b(l) + d(0)). These attempted events take time 
At = 1/n on average, therefore time is increased by At. 
Variation of the background can be handled in the same 
way within a smaller (I x I) region whose sites affect the 
CP at the given time. Namely, the values of g(r) are 
redrawn at randomly chosen Atfl 2 sites. 

The average values defined above [see Eqs. are 
determined for discrete time values chosen equidistantly 
in a logarithmic scale. 

IV. RESULTS FOR FROZEN BACKGROUND 

First we concentrate on a system where the random 
background is frozen, i.e. / = 0. Varying the propor- 
tion of good sites (P), we have made a systematic MC 
analysis for fixed death and birth rates. For the present 
parameters (see Fig. |l|) the good sites form isolated clus- 
ters [having P < P pcrc = 0.5926 (percolation threshold)] 
and infinite spreading is possible across the bad sites. 
The results summarized in the subsequent figures agree 
quantitatively with those found by Moreira and Dickman 
H on a diluted lattice for P > P per c- 




time [MCS] 

FIG. 1. Survival probabilities as functions of time at dif- 
ferent P values (indicated by labels) for b(l) = 1, 6(0) = 0.5, 
d(l) = 0.25, and d(0) = 0.5 . 

Figure [l] shows some typical time-dependence functions 
of survival in a log- log plot. The statistical errors are 
comparable to the line width due to the large number of 
trials (e.g. M — I0 8 for P — 0.22) . In the active region 
(P > P c ), S(t) tends to a constant value as expected p4[ . 
In the subcritical region, S(t) can be well approximated 
by a power law. Unfortunately, convergency toward the 
limit value becomes extremely slow in the vicinity of the 
critical point, and this prevents the accurate determi- 
nation of P c . Division between the active and inactive 
regions becomes more visible when considering the func- 
tion N(t) (see Fig. §). 
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FIG. 2. Average number of individuals vs. time for the 
same simulations as in Fig. |l|. 

In the active region, the surviving individuals occupy a 
compact patch whose radius increases linearly with time. 
Consequently, N(t) becomes proportional to t 2 for suffi- 
ciently long times as indicated in the Fig. ^. By contrast, 
when P < P c , then N(t) tends to a power law decay with 
an exponent depending on P. 
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FIG. 3. Log-log plot of mean-square distance of individuals 
from the origin as a function of time. Parameter values are 
the same as in Fig. jjj. 

Figure |^ shows the time evolution of the mean square 
distance of surviving individuals as defined by Eq. (||). 
In agrement with the expectations R 2 (t) oc t 2 for suffi- 
ciently long times if P > P c . In the subcritical region, 
the increase of R 2 (t) becomes significantly slower than 
that found for the homogeneous system. In the homoge- 
neous system, the surviving individuals perform random 
walks independently of each other, therefore R 2 (t) oc t in 
the subcritical region. In the present case, however, our 
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data indicate a power law increase: R 2 (t) oc t x , where 
the exponent A(P) < 1. This indicates localization of 



spreading as, it will be explained in Section VI. 

Evaluation of N(t) involves those (stoped) trials which 
have reached the empty state before the given time t. 
Consequently, the ratio N(t)/S(t) expresses the aver- 
age number of individuals in the surviving trials. In the 
subcritical region, the time-dependence of this quantity 
tends to logarithmic increase as shown in a lin-log plot 
(see Fig. |). 
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FIG. 4. Time dependence of the average number of indi- 
viduals in the surviving trials at P — 0.22, 0.23, and 0.235. 

In the classical theory of percolation (reviewed in Refs. 
jl8],[l9]]), a geometrical feature of the clusters of size s 
is characterized by an average cluster radius i? per c( s )- 
Below the percolation threshold R pC rc(s) oc s p for the 
large s limit. According to MC simulations p8| , p0| , ^l| , 
p ~ 2/3, while the "self-avoiding walk" prediction fffi] 
gives p ~ 3/4. In general, p > 1/2 values are characteris- 
tic for ramified clusters . On the contrary, above the 
percolation threshold the finite size clusters are compact, 
that yields p — 1/2. 

Following this analogy, we examine the relation be- 
tween N(t)/S(t), the actual size of population in the 
surviving trials, and the average radius of the area they 
occupy at a given time in the subcritical region. Figure ^ 
shows this relation on a log- log plot. In the close vicinity 
of the critical point (at P = 0.239) we find power law 
behavior with an exponent 1.66(6). This behavior indi- 
cates that the species spreads over a ramified cluster in 
the whole time interval we could study the system. 

In the subcritical region, however, we can observe two 
different regions. In the transient (short-time) region, the 
increase of R 2 {t) is faster than at the critical point. This 
indicates that the offspring move away from the origin 
along preferred pathes formed by good sites. Conversely, 
the long-time behavior is dominated by those trials where 
the individuals stay on a compact patch. 
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FIG. 5. R 2 (t) vs. N(t)/S(t) for P = 0.22 (dashed-dotted), 
0.23 (dashed), 0.235 (solid), and 0.239 (dotted line). 

In the next section we show that these features funda- 
mentally change when we allow the background change 
over time. 



V. RESULTS FOR TIME-DEPENDENT 
BACKGROUND 

Applying time-dependent backgrounds, we determined 
the same quantities [S(t), N(t), and R 2 (t)] as before. 
We use the same birth and death rates when varying the 
value of P for a fixed rate (/) of background change. 
Our results estimate the critical value for the ratio of 
good sites at P c = 0.2680(1). 
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FIG. 6. Log-log plot of survival probabilities versus time 
at different P values (as indicated by labels) for 6(1) = 1, 
6(0) = 0.5, d(l) = 0.25, d(0) = 0.5, and / = 0.01 . 

Figure |^ shows that the survival probability S(t) van- 
ishes exponentially when P < P c . In the active region 
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(P > P c ), however, S(t) tends to a constant value. No- 
tice, that the convergence is significantly faster here than 
at the frozen background (see Fig. |l|). 

Data for P = 0.268 represent the behavior of the sys- 
tem at the critical point. A numerical analysis confirms 
that the survival probabilities can be well approximated 
by a power law, S(t) cx t~ s for sufficiently long times. 
Numerical fitting gives S = 0.45(2) in agreement with 
the exponent found in DP transitions in homogeneous 
systems [|,|. 

When considering N(t), we can also distinguish be- 
tween three different behaviors (see Fig. [?]). Above the 
critical point (P > P c ), N(t) tends towards quadratic 
time dependence. Below the critical point, N(t) de- 
creases exponentially. At the critical point, however, our 
MC data indicate power law behavior, i.e., N(t) cx f 
with r] — 0.23(1). This numerical value of r\ agrees with 
the characteristic exponent of DP transitions [||J|] . 




time [MCS] 

FIG. 7. Average number of individuals vs. time for the 
same simulations as in Fig. |^. 

In summary, all the above-mentiend numerical results 
support that contact process exhibits the same features 
in time-dependent random environments as in homoge- 
neous background. Similarly to S(t) and N(t), R 2 (t) also 
supports this conclusion (data not shown). 

These results have been gained at a fixed rate of 
changes (/ = 0.01). Preliminary results suggest, how- 
ever, that the same conclusions holds for other values 
of / as well. General features of spreading remain un- 
changed due to the robustness of DP transition. 

Here it is worth mentioning that we have also stud- 
ied the concentration of individuals (as well as its fluc- 
tuation) in the stationary state for some parameter val- 
ues. These analyses support that the variation of both 
quantities is consistent with the expectations, e.L, the 
species concentration and its fluctuation exhibit power 
law behavior with exponents characteristic for the two- 
dimensional DP systems ]T^ , p^ ,^,pt . 



VI. CONCLUSIONS AND SUMMARY 

We have studied a generalized version of contact pro- 
cess on a square lattice. The background in the present 
model provides a continuous transition from a homoge- 
neous system to a randomly fluctuating environment, 
including (frozen) diluted lattices. The environment is 
inhomogeneous, consisting of two types of sites (good 
and bad) that affect the rates of birth and death of the 
species. In this model, spreading from a cluster of good 
sites to another one is permitted across the bad sites 
and/or bridges (of good sites) created occasionally by 
background fluctuations. 

Our analyses have been rescricted to the considera- 
tion of some typical features, because the simulations are 
rather time-consuming. The results confirm the theo- 
retical expectations as well as the general picture drawn 
by previous authors who considered CP on diluted 
lattices. 

Obviously, an inhomogeneous environment contains 
areas that provide better (or worse) conditions for the 
survival relative to the average quality of the habitat. 
This variation becomes particularly important in the 
vicinity of the transition point that separates the active 
and inactive regions for long-time behavior. In the sub- 
critical region, the surviving individuals are constrained 
to remain within isolated patches. Conversely, in the ac- 
tive region, favorable areas can form a percolating clus- 
ter, that has infinite extension, and sustains infinite sur- 
vival. 

To explain the main results, we briefly remind the 
reader of a simple calculation suggested by Noest J7| . In 
this description wc characterize the favorable patches by 
size s and assume that their probability decreases expo- 
nentially, i.e., n s cx e~ As . Furthermore, it is also assumed 
that the average over many trials (or patches of size s) 
yields an exponentially decreasing number of individuals 
over time, that is, m s (t) cx e _ *' Ts where the average sur- 
vival time can be approximated as r s cx e Bs if initially all 
sites are occupied. Noest has shown that long-time be- 
havior can be well approximated in this case by a power 
law decrease, that is, m(i) = ^ s sn s m s (t) cx t~ A l B . 
This prediction is derived from a maximum likelihood 
estimation that implies the existence of a typical cluster 
size giving the dominant contribution to m(t). According 
to this approach, the size of dominant clusters increases 
logarithmically with time t. 

The adoption of this approach is not straightforward 
to the present situation where we consider the spreading 
from a single seed. In our case, the species can die out 
before achieving homogeneous distribution, while homo- 
geneous initial state is assumed in the former calculation. 
We can assume, however, that the species spreads over 
the whole cluster or dies out within a transient time, and 
the probabilities of these outcomes are independent of 
cluster size. If this hypothesis holds, than the approach 
of Noest can be applied for the analysis of long-time be- 
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and the shape-dependence of r s . 

Numerical analysis of the extinction process is difficult 
(in a quenched environment) because decrease of S(t) 
and N(t) slowes down in the vicinity of the critical point. 
In practice this implies persistence over long time peri- 
ods before extinction. The slow extinction process, how- 
ever, requires rather specific environmental conditions. 
Even a small degree of temporal fluctuation results in 
an extinction process that is analogous to DP transition 
on large scales. Our mean-field calculations and pre- 
liminary simulations have indicated that the transition 
point is strongly affected by the rate / of background 
change. Further analyses are required to quantify this 
phenomenon for different birth and death rates. 
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Some predictions of this rough method have been con- 
firmed by our MC simulations. Namely, S(t) and N(i) [10 
exhibit a dominant power law decrease in the subcritical 
region, while N(t)/S(t) (in proportion to the dominant [H 
cluster size) increases logarithmically. Furthermore, in 
the active region, the algebrically decaying contribution 
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of S(t) comes from extinction on (compact) finite size 
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clusters whose probability distribution satisfies the above 
assumption 21 h|. 

Noest 's approach is focused on compact clusters where 
s tx R 2 . Consequently, N(t)/S(t) oc R 2 (t) is expected 



when the system behavior is dominated by a typical clus- 
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ter size. This behavior can be observed for long times far rj« 
below the critical point. Figure |^ indicates, however, 

that these predictions (and assumptions) fail in the close r^gi 

vicinity of the critical point. In other words, the ramified j^g] 

clusters give relevant contribution to S(t) and N(t) dur- pg 
ing a transient time that increases when approaching the 

critical point. In this case, we need a more sophisticated [21] 

description allowing the emergence of ramified clusters 22 



G 



